# konrad.smolinski@gmail.com
# date: 		26/11/2010
# last update: 	12/05/2011
#
# info:
#	this file works with an illustration from rcprobit4.pdf and rcprobit1.nb

#	(adc) mail note: 13/05/2011
#	In my initial grid experiments I have used a binary instrument 
#	z in {-1,1}  
#	a1=0, a2=1, b1=1, b2=1  and d1=1 (not 2).  
# ------------------------------------------------------------------
# Call for libraries:
rm(list=ls())
source("./rcLibSrc.R")
# ---- Set Up
# following 'rcprobit4.pdf'
	K <- 2 

# PARAMETERS
	a <- c(0,1)
	b <- c(0,1)
	d <- c(1,3^(-1/2),3^(-1/2),3^(-1/2))
	zval <- c(-0.5,0.5)
	R <- length(zval)
	xthres <- c(0)
	xval <- c(0,1)

# integration accuracy
	acc <- 40

# Validate :
	dgp <- rcDgp(a,b,d,xthres,xval,zval)
	dgp

# SIMPLE EXAMPLE 
# this part is done once:
	probIdx <- rcLeaveOutToProbs(K)
	probMat <- matrix(0,ncol=R,nrow=length(probIdx))

	for(r in 1:R){
		for(i in 1:length(probIdx)){
			 probMat[i,r] <- sum(dgp[[r]][probIdx[[i]]]) 
		}
	}
	rhs <- apply(probMat,1,max)


# integrand over Base Components - this part depends on initial parameters
	loIdx <- rcLeaveOutIdx(K)
	lhs <- rep(0,length(loIdx))

# parameters enter here:
	intVec <- rcBaseIntegrate(a=c(0,1.05),b=c(0,1),xval=xval)
	for(i in 1:length(loIdx))	lhs[i] <- sum(intVec[loIdx[[i]]])
	lhs

lhs > rhs

all(lhs > rhs)
# END OF TEST
# ------------------------------------------------------------------------------
# GRID
	a0seq <- seq(-1.6,1.6,len=321)
	a1seq <- seq(-1.0,2.6,len=321)
	aGrid <- expand.grid(a0seq,a1seq)

	b1 <- 0
	b2 <- 1

	res <- rep(FALSE,nrow=nrow(aGrid))
	resInqMat <- matrix(FALSE,nrow=nrow(aGrid),ncol=length(lhs))
	qs <- smolyak.quad(2,40)

l <- 1
for(g in 1:nrow(aGrid)){
	lhs <- rep(0,length(loIdx))
	intVec <- rcBaseIntegrate(a=c(aGrid[g,][[1]],aGrid[g,][[2]]),b=c(b1,b2),xval=xval,quadPt=qs)

	for(i in 1:length(loIdx))	lhs[i] <- sum(intVec[loIdx[[i]]])

	resInqMat[g,] <- lhs - rhs
 	res[g] <- all(resInqMat[g,] > 0)
	cat(".")
}
	
	set <- aGrid[res,]



save(set,dgp,a0seq,a1seq,a,b1,b2,d,zval,xval,xthres,file='../Data/rcEx5-1-Data.RData')
load(file='../Data/rcEx5-1-Data.RData')
# plot
polyAlpha <- list()
sv <- T
if(sv) pdf("../Graphs/rcEx5-1.pdf",width=10,height=6)
#postscript("../Graphs/rxEx1-1.eps",width=10,height=6)

	par(mfrow=c(1,1))

	polyAlpha <- ahull(set,alpha=0.1)
	plot(polyAlpha,main=paste("b[1]=0","b[2]=1",sep="")
	,xlab=expression(a[1])
	,ylab=expression(a[2])
	,wlines='vor',col=c(1,2,3,3,0,2))
	
	abline(v=0,col="red",lty=1)
	abline(h=1,col="red",lty=1)
	points(0,1,pch=19,cex=1.2,col=2)


if(sv) dev.off()